# 灰色模型预测

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

def gm(dirt,stname,num,title,s):
    """dirt：数据文件绝对地址。比如，C:\\Users\\sddbook\\Desktop\\gmtest.xlsx。
    stname：Excel中哪个表，比如Sheet2。
    num：数据总个数。title：数据所在列的字段名，便于计算机确定是哪一列数据.s:预测个数"""
    data = pd.read_excel(dirt, sheetname=stname)
    data = np.array(data[title])
    lens = len(data)  # 数据量

    # 数据检验
    ## 计算级比
    lambds = []
    for i in range(1, lens):
        lambds.append(data[i-1]/data[i])
    ## 计算区间
    X_min = np.e**(-2/(lens+1))
    X_max = np.e**(2/(lens+1))
    ## 检验
    is_ok = True
    for lambd in lambds:
        if (lambd < X_min or lambd > X_max):
            is_ok = False
    if (is_ok == False):
        print('该数据未通过检验')
    else:
        print('该数据通过检验')

    # 构建灰色模型GM(1,1)
    ## 累加数列
    data_1 = []
    data_1.append(data[0])
    for i in range(1, lens):
        data_1.append(data_1[i-1]+data[i])
    ## 灰导数及临值生成数列
    ds = []
    zs = []
    for i in range(1, lens):
        ds.append(data[i])
        zs.append(-1/2*(data_1[i-1]+data_1[i]))
    ## 求a、b
    B = np.array(zs).reshape(lens-1,1)
    one = np.ones(lens-1)
    B = np.c_[B, one]  # 加上一列1
    Y = np.array(ds).reshape(lens-1,1)
    a, b = np.dot(np.dot(np.linalg.inv(np.dot(B.T, B)), B.T), Y)
    print('a='+str(a))
    print('b='+str(b))

    # 预测
    def func(k):
        c = b/a
        return (data[0]-c)*(np.e**(-a*k))+c
    data_1_hat = []  # 累加预测值
    data_0_hat = []  # 原始预测值
    data_1_hat.append(func(0))
    data_0_hat.append(data_1_hat[0])
    for i in range(1, lens+s):  # 多预测s次
        data_1_hat.append(func(i))
        data_0_hat.append(data_1_hat[i]-data_1_hat[i-1])
    print('预测值为：')
    for i in data_0_hat:
        print(i)

    # 模型检验
    ## 预测结果方差
    data_h = np.array(data_0_hat[0:num]).T
    sum_h = data_h.sum()
    mean_h = sum_h/lens
    S1 = np.sum((data_h-mean_h)**2)/lens
    ## 残差方差
    e = data - data_h
    sum_e = e.sum()
    mean_e = sum_e/lens
    S2 = np.sum((e-mean_e)**2)/lens
    ## 后验差比
    C = S2/S1
    ## 结果
    if (C <= 0.35):
        print('1级，效果好')
    elif (C <= 0.5 and C >= 0.35):
        print('2级，效果合格')
    elif (C <= 0.65 and C >= 0.5):
        print('3级，效果勉强')
    else:
        print('4级，效果不合格')

    #画图
    # plt.figure(figsize=(9, 4), dpi=100)
    # x1 = np.linspace(1, 7, 7)
    # x2 = np.linspace(1, 12, 12)
    # plt.subplot(121)
    # plt.title('x^0')
    # plt.plot(x2, data_0_hat, 'r--', marker='*')
    # plt.scatter(x1, data, marker='^')
    # plt.subplot(122)
    # plt.title('x^1')
    # plt.plot(x2, data_1_hat, 'r--', marker='*')
    # plt.scatter(x1, data_1, marker='^')
    # plt.show()
#gm('C:\\Users\\sddbook\\Desktop\\gmtest.xlsx','Sheet3',26,'test'，5)